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Quasi-Periodic Oscillations in Relativistic Tori 
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Motivated by recent interesting work on p-mode oscillations in axisymmetric hydrodynamic black- 
hole tori by Rezzolla, Zanotti, and collaborators, I explore the robustness of these oscillations by 
means of two and three-dimensional relativistic hydrodynamic and MHD simulations. The primary 
purpose of this investigation is to determine how the amplitudes of these oscillations are affected by 
the presence of known instabilities of black-hole tori, including the Papaloizou-Pringle instability 
(PPI) and the magneto-rotational instability (MRI). Both instabilities drive accretion at rates above 
those considered in Rezzolla’s work. The increased accretion can allow wave energy to leak out of 
the torus into the hole. Furthermore, with the MRI, the presence of turbulence, which is absent 
in the hydrodynamic simulations, can lead to turbulent damping (or excitation) of modes. The 
current numerical results are preliminary, but suggest that the PPI and MRI both significantly 
damp acoustic oscillations in tori. 


I. INTRODUCTION 


“Diskoseismology” jl] is a relatively new subfield of 
accretion disk physics. Often diskoseismic oscillations 
have been overlooked as uninteresting because, to be 
observable from the outside, they must be trapped in 
a particular region of the disk or have a global pattern. 
For a Keplerian disk in a Newtonian potential, there 
are no global modes and relatively few trapped modes. 
Relativistic gravity helps the situation somewhat by 
providing a richer set of trapped modes. Nevertheless, 
even in general relativity, Keplerian thin disks have 
only a few trapped modes, which mostly occur in very 
restricted regions of the disk. 

Non-Keplerian disks, on the other hand, have a 
much richer set of diskoseismic modes available for 
two fundamental reasons. First, in contrast to Keple¬ 
rian disks, which in principle extend to infinity, non- 
Keplerian disks can be constructed to have a finite size 
and can thus act as a single resonant cavity. Second, 
unlike Keplerian disks, non-Keplerian disks necessar¬ 
ily have radial pressure gradients that can provide an 
additional restoring force (along with the centrifugal 
force and gravity) in support of oscillations. 

In this work I focus on inertial-acoustic (p-mode) os¬ 
cillations in finite tori. Two interesting consequences 
of such oscillations have previously been identified. 
First, they can cause periodic changes in the mass 
quadrupole moment of the disks. For massive tori 
of nearly nuclear densities, this may produce gravi¬ 
tational radiation detectable by LIGO, especially for 
sources located within our galaxy |2]. Such massive 
tori could form from the gravitational collapse of a 
massive rotating star or as an intermediate state in 
a binary neutron star merger. These oscillations may 
also be important in explaining quasi-periodic oscil¬ 
lations (QPOs), particularly the harmonically related 
high-frequency QPOs (HFQPOs) seen in some black- 
hole candidate low-mass X-ray binaries (LMXBs) flj. 

My goal here is to extend previous numerical sim¬ 
ulations of oscillating tori II] to include magnetic 


fields and nonaxisymmetric perturbations. I proceed 
in m by reviewing the equations of general relativis¬ 
tic magnetohydrodynamics as solved in the numerical 
code used in this work. The code itself is briefly de¬ 
scribed in m In M\ I present the results of this 
study. I conclude with some final thoughts in El 


II. GENERAL RELATIVISTIC MHD 
EQUATIONS 


This work uses a form of the general relativistic 
MHD equations similar to De Villiers and Hawley 
(B). Nevertheless, it is worth explicitly writing these 
equations out for clarity. In flux-conserving form, the 
conservation equations for mass, internal energy, and 
momentum and the induction equation for magnetic 
fields are: 
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where g is the determinant of the 4-metric, W = 
yf^gu 1 is the relativistic boost factor, D = Wp is 
the generalized fluid density, P = (P — 1 )E/W is the 
fluid pressure, V 1 = u % /uf is the transport velocity, 
Si = Wphui is the covariant momentum density, and 
E = We = Wpe is the generalized internal energy 
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density. There are two representations of the mag¬ 
netic field in these equations: B M is the magnetic field 
4-vector (||B|| 2 = g^B^B’') and 

B = W{B^ - B a V») (5) 

is the divergence-free (dB l /dx l = 0), spatial (B° = 0) 
representation of the field. The time component of the 
magnetic field B° is recovered from the orthogonality 
condition B^u^ = 0. 

These equations are evolved in a Kerr-Schild polar 
coordinate system (t, r , 0, 4>), although all of the sim¬ 
ulations in this work assume a non-rotating (a = 0) 
Schwarzschild black hole. The computational advan¬ 
tages of the “horizon-adapted” Kerr-Schild form of 
the Kerr metric are described in Papadopoulos and 
Font [§] and Font et al. Q. The primary advantage is 
that, unlike Boyer-Lindquist coordinates, there are no 
singularities in the metric terms at the event horizon. 
This is particularly important for numerical calcula¬ 
tions as it allows one to place the grid boundaries in¬ 
side the horizon, thus ensuring that they are causally 
disconnected from the rest of the flow. 


III. NUMERICAL METHOD 

These simulations are carried out using the numeri¬ 
cal code Cosmos!—|-, a massively parallel, multidimen¬ 
sional (one, two, or three dimensions), adaptive-mesh, 
magnetolrydrodynamic code for evolving both New¬ 
tonian and relativistic flows. Cosmos++ employs a 
time-explicit, operator-split, finite volume discretiza¬ 
tion method with second-order spatial accuracy. A 
detailed description of Cosmos-1—b, including test re¬ 
sults, will be presented in an upcoming paper. 

The simulations are initialized with the analytic so¬ 
lution for an axisymmetric torus with constant spe¬ 
cific angular momentum l. For the initialization, we 
assume an isentropic equation of state P = np T , 
although during the evolution, the adiabatic form 
P = (F— 1)E/W is used to recover the pressure when 
solving equations @ and ©. We set V = 4/3 and 
ft = 0.0229 (in G = c = 1 units). 

This work presents both two-dimensional (axisym¬ 
metric) and three-dimensional simulations. The two- 
dimensional simulations are carried out on a grid ex¬ 
tending from 0.98 vbh < r < r max and 0 < 0 < 7r, 
where xbh = 2GM/c 2 is the radius of the black-hole 
horizon. The three-dimensional simulations include 
the full azimuthal range 0 < <fi < 2i r. We choose 
r m ax to be about twice the initial outer radius r ou t of 
the torus. In order to increase the resolution inside 
the torus, we replace the radial coordinate r with a 
logarithmic coordinate x\ = 1 + ln(r/rs#) and the 
angular coordinate 0 with a coordinate X 2 satisfying 
6 = X2+\sm(2x2). The grid is resolved with 128x 128 
zones in 2D, giving a radial spacing of 0.04rG near the 


horizon and 0 .6ro near the outer boundary. In the 
present work, the 3D simulations are limited to 64 3 
zones. 

In the “background” regions not determined by the 
initial torus solution, we initialize the gas following 
the spherical Bondi accretion solution Jgj. We fix the 
parameters of this solution such that the rest mass 
present in the background is negligible compared to 
the mass in the torus. The outer radial boundary 
is held fixed with the analytic Bondi solution for all 
evolved fields. The inner radial boundary uses out¬ 
flow ( V r < 0) boundary conditions. Data are shared 
appropriately across angular boundaries. 


IV. RESULTS 

Results are presented in geometrized units (G = c = 
1) with units of length parameterized in terms of the 
gravitational radius of the black hole, re = GM/c 2 . 


A. 2D Axisymmetric Hydro 

It is instructive to begin this work by reproduc¬ 
ing the results of a previously published axisymmet¬ 
ric hydrodynamic simulation of an oscillating torus. 
This will be useful for illustrating how p-mode oscil¬ 
lations are generally manifested in axisymmetric hy¬ 
drodynamic tori and facilitate easier comparison with 
the subsequent MHD and non-axisymmetric results. 
This is also important as a validation of the new code 
used here. I choose Model (a) from Zanotti et al. 
[2] for this purpose. In this model, the specific an¬ 
gular momentum, which is constant throughout the 
torus, is l/M = 3.8. The torus just fills the largest 
closed equipotential surface so that the inner torus 
boundary corresponds with the location of the cusp 
r m = r cusp = 4.576. The outer radius for such a torus 
is r out = 15.889. The density center of this torus is 
located at r cen ter = 8.352; the orbital period at this 
radius is t or t, = 151.7. This will serve as a reference 
dynamical timescale for all these simulations. To ex¬ 
cite the resonant acoustic mode I use an initial radial 
velocity kick. For convenience this velocity is set at 
a small fraction (rj = 0.12) of the background Bondi 
inflow velocity. 

I use the L 2 norm of the rest mass density, defined 
as I HI 2 = Piji to characterized the global 

oscillatory behavior of the torus. The time history and 
Fourier power spectrum of this quantity are shown in 
Figures n and [2 respectively. The power spectrum, in 
particular, reveals the rich harmonic structure char¬ 
acteristic of diskoseismic modes ||9]. Rezzolla et al. 
also showed that this structure is consistent with pre¬ 
dictions of linear perturbation analysis of vertically 
integrated relativistic tori. 
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FIG. 1: Time history of the L 2 norm of the rest-mass 
density. The vertical scale has been adjusted to the initial 
value of ||p|| 2 . 



FIG. 2: Power spectra of the L 2 norm of the rest-mass 
density. The spectra have been arbitrarily normalized 
such that the strongest frequency bin has a value of 
100. The frequency is set assuming a black hole mass 
Mbh = 2.5Mq. 


Although the p-mode of interest in this work is an 
acoustic wave within the disk, the torus parameters for 
this model are such that this wave periodically pushes 
material over the cusp and out of the torus. This ma¬ 
terial is accreted into the black hole on a dynamical 
timescale. Therefore, one manifestation of the p-mode 
in (marginally stable) tori is a periodic fluctuation in 
the black-hole mass accretion rate as shown in Fig¬ 
ure El The Fourier power spectrum of the accretion 
history (Figure EJ reveals the same fundamental fre¬ 
quency and first overtone as Figure El confirming the 
connection. 



FIG. 3: Rest mass accretion rate normalized by the initial 
mass of the torus. 



FIG. 4: Power spectra of the mass accretion rate. The 
spectra have been arbitrarily normalized such that the 
strongest frequency bin has a value of 100. The frequency 
is set assuming a black hole mass Mbh = 2.5 Mg. 


B. 3D Hydro 

I now drop the restriction of axisymmetry and con¬ 
sider a three-dimensional hydrodynamic torus. The 
concern here is that such a structure is known to be 
unstable to low-order non-axisymmetric modes (the 
Papaloizou-Pringle instability or PPI). Note that this 
simulation uses a slightly different initial setup than 
the previous one. First, the velocity perturbation is 
slightly smaller in this case (77 = 0.06 instead of 0 . 12 ). 
In my testing I have found that the amplitude of this 
perturbation has little effect on the results, so this 
difference shouldn’t be important. This simulation 
also uses a different set of parameters for the back¬ 
ground Bondi solution. These differences account for 
the much larger amplitude of oscillations in the ac¬ 
cretion rate for the three-dimensional model (Figure 
El) Nevertheless, the frequency of the oscillations in 
both cases is the same, indicating that the different 
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FIG. 5: PPI mode growth. 


background treatments do not effect the behavior of 
the p-rnode oscillations. 

It is useful for this simulation to add a diagnostic 
to track the growth of the main PPI modes. We ex¬ 
tract the m = 1 and 2 modes by computing azimuthal 
density averages as 10 ] 

Re[/c m (r)] = / p(r, 7 r/ 2 , </>) cos(ra</>)d(/> , ( 6 ) 

Jo 

p2n 

lm[k m (r)]= / p(r, 7 r/ 2 , 0 )sin(m 0 )d 0 . (7) 

Jo 

The power in mode m is then 



- 3.0 - 2.0 - 1.0 0.0 1.0 2.0 3.0 


FIG. 6: Density in the midplane of the torus. This plot is 
shown using the logarithmic radial coordinate. 

in the internal structure of the torus. The damping 
is probably a result of leaking wave energy into the 
black hole along the accretion stream. 

C. 2D Axisymmetric MHD 


fm = 


Tmax Tmin J rmi r 

+ {In# m (r)]} 2 ) 


In 


({Re[fc m (r )]} 2 ( 8 ) 


dr , 


where r m i n and r max are the approximate inner and 
outer edges of the disk. In Figure 0 we show the 
growth of these two modes as a function of time. As 
expected, there is significant mode growth for this 
slender torus. Furthermore, this mode growth does 
not appear to saturate before the end of the simu¬ 
lation. Figure El shows the midplane density of this 
model at t/t or b = 33. An m = 1 “planet” and inflow¬ 
ing spiral wave are both apparent. 

The question of interest now is what effect this PPI 
growth has on the p-mode oscillations. From Figure 
131 it is clear that the PPI growth coincides with a sig¬ 
nificant increase in the mass accretion rate, due to 
the angular momentum transport of the spiral waves 
generated in the torus. Despite the very significant 
changes to the structure of the torus, the fundamen¬ 
tal p-mode is able to survive at least to the end of this 
simulation, as is apparent in Figure 0 However, the 
power in this mode is significantly reduced and it is 
smeared out in frequency space as shown in Figured 
The frequency broadening is probably due to changes 


I now explore the effect of adding initially weak 
poloidal magnetic field loops to the axisymmetric 
torus investigated in The presence of the 

poloidal field triggers the violent growth of the 
magneto-rotational instability [MRl IllLlT^] . The MRI 
facilitates angular momentum transport through mag¬ 
netohydrodynamic turbulence. Similar to - II V III the 
focus here is to determine whether the p-rnode oscil¬ 
lations seen in the 2D hydrodynamic simulation can 
survive in the presence of a known instability, in this 
case the MRI. 

The initial magnetic field vector potential is @| 

_ f Pmin ) fol" P fl. Pmin 5 /g\ 

0 ~ \ o for p < p min . 1 ’ 

The non-zero spatial magnetic field components are 
then B r = — dgA$ and B e = d r A These poloidal 
field loops coincide with the isodensity contours of the 
torus. The constant k is normalized such that initially 
/3 = P/(\\B\\ 2 /8tt) > 100 throughout the torus. The 
parameter = 0.5 *pcenter is used to keep the field 
a suitable distance inside the surface of the torus. 

The presence of MRI turbulence dramatically in¬ 
creases the mass accretion rate compared to the hy¬ 
drodynamic torus (see Figurc|«3l and also redistributes 
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FIG. 7: Logarithm of density for the initial and final torus 
states for the 2D MHD simulation. 

the disk material into a broader, more radially ex¬ 
tended structure as shown in Figure0 It is important 
to note that axisymmetric MRI simulations are sus¬ 
ceptible to a particularly violent form of the poloidal 
MRI, called the channel solution [llj, which can itself 
lead to highly episodic mass accretion [cf. Hi , appar¬ 
ent in Figure 0 

Comparisons of the power spectra of the L 2 density 
norms (Figure^) and mass accretion rates (Figure0 
suggest that the violent nature of the channel solution 
strongly damps the p-mode oscillations. The accretion 
spectrum of this simulation shows no power in the 
fundamental mode or its overtones. The L 2 density 
norm spectrum shows some power in the fundamental 
mode, although the amplitude of the oscillation decays 
with time, as apparent in Figure 0 It remains to be 


seen if this strong damping is restricted to the onset of 
the MRI or is a more general feature of these turbulent 
tori. 

V. SUMMARY 

An obvious extension of this work will be to perform 
3D MHD simulations. There are two potentially im¬ 
portant differences between magnetized axisymmetric 
(2D) and nonaxisymmetric (3D) tori. One is the fact 
that the very violent channel solution of MRI (seen 
in Figure 0 is itself susceptible to a nonaxisymmet¬ 
ric instability that destroys its coherence [lfj. As a 
consequence, accretion from a nonaxisymmetric torus 
is much less episodic than the equivalent 2D torus 
0 The second is that the MRI is unable to sustain 
itself through dynamo action in axisymmetry. The 
sustained activity of the MRI in three-dimensions will 
be important in assessing the importance of p-mode 
oscillations in magnetized relativistic tori. 

If acoustic modes are indeed strongly damped in 
turbulent MRI disks, then the puzzle of QPOs in 
black-hole X-ray binaries becomes more problematic. 
The MRI only requires that a disk be differentially 
rotating at a rate that decreases with radial distance 
and that it have a weak magnetic field. Both of these 
conditions are thought to be generally met in astro- 
physical accretion disks. In light of this and the re¬ 
sults presented here, it appears that p-mode oscilla¬ 
tions may not be a viable mechanism for generating 
observable QPOs in realistic accretion disks, but 3D 
MHD simulations will be necessary to confirm this. 

Animated movies of these results are available at 
http://www.physics.ucsb. edu/~ fragile/research. html. 
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